#ifndef _DIPOLARITY_HH_
#define _DIPOLARITY_HH_

#include "fastjet/PseudoJet.hh"
#include <vector>
#include <cstdlib>

class Dipolarity
{
  // signed phi distance function
  double signedphidist(double phia, double phib)
  {
    double deltaphi=std::abs(phia-phib);
    double phiL=phia, phiS=phib;
    if (phib>phia) { phiL=phib; phiS=phia; }
    
    if (6.28318531-phiL+phiS<deltaphi)
      {
	deltaphi=6.28318531-phiL+phiS;
	if (phia<phib) deltaphi*=-1;
	return deltaphi;
      }
    
    if (phia>phib) deltaphi*=-1;
    
    return deltaphi;
  }


  double operator()(const fastjet::PseudoJet& subjet1, 
		    const fastjet::PseudoJet& subjet2, 
		    const std::vector<fastjet::PseudoJet>& subjets)
  {
    double y, phi, dist;
    double dipolarity=0.0, u, yint, phiint;
    double y1 = subjet1.rap(), y2 = subjet2.rap();
    double phi1 = subjet1.phi(), phi2 = subjet2.phi();
    double radius = std::sqrt(subjet1.plain_distance(subjet2));
    fastjet::PseudoJet totjet(0,0,0,0);
    
    for (int k=0; k<(int)subjets.size(); k++)
      {
	totjet+=subjets[k];
	y = subjets[k].rap();
	phi = subjets[k].phi();
	u = (y-y1)*(y2-y1)+signedphidist(phi1,phi)*signedphidist(phi1,phi2);
	u /= ((y2-y1)*(y2-y1)+std::pow(signedphidist(phi1,phi2),2));
	if (u>1.0) u=1.0;   if (u<0.0) u=0.0;
	yint = y1 + u*(y2-y1);
	phiint = phi1 + u*signedphidist(phi1,phi2);
	if (phiint>6.28318531)   phiint-=6.28318531;
	if (phiint<0.0)          phiint+=6.28318531;
	dist = std::sqrt(std::pow(signedphidist(phiint,phi),2)+std::pow(yint-y,2));
	dipolarity+=subjets[k].perp()*dist*dist;
      }
    
    dipolarity/=totjet.perp();
    dipolarity/=(radius*radius);
    return dipolarity;
  }
  
  
};
#endif
